#include <iomanip>

#include "snp_call.hpp"
#include "util.hpp"
#include "option_parse.hpp"
#include "call_rates.hpp"

void multiply_seq_snp_calls_by_priors(prog_options &options, seq_call_vector_t &snps)
{
    for(unsigned int i = 1; i < snps.size(); i++)
    {
        SeqSnpCall &snp = snps[i];
        
        if(!options.diploid)
            for(unsigned int j = 0; j < snp.probs.size(); j++)
                snp.probs[j] *= options.haploid_priors[cswap[snp.seq[0]]][cswap[snp.blocks[j].base]];
        else
            for(unsigned int j = 0; j < snp.probs.size(); j++)
                snp.probs[j] *= options.diploid_priors[cswap[snp.seq[0]]][cswap[snp.blocks[j].base]];

        // normalize score.
        double sum = 0.0;
        for(unsigned int j = 0; j < snp.probs.size(); j++)
            sum += snp.probs[j];
        for(unsigned int j = 0; j < snp.probs.size(); j++)
            snp.probs[j] /= sum;
    }
}

void multiply_col_snp_calls_by_priors(prog_options &options, color_call_vector_t &snps)
{
    for(unsigned int i = 1; i < snps.size(); i++)
    {
        ColorSnpCall &snp = snps[i];


        if(!options.diploid)
            for(unsigned int j = 0; j < snp.probs.size(); j++)
                for(unsigned int k = 0; k < snp.blocks[j].size(); k++)
                    snp.probs[j] *= options.haploid_priors[cswap[snp.seq[k]]][cswap[snp.blocks[j][k].base]];
        else
            for(unsigned int j = 0; j < snp.probs.size(); j++)
                for(unsigned int k = 0; k < snp.blocks[j].size(); k++)
                    snp.probs[j] *= options.diploid_priors[cswap[snp.seq[k]]][cswap2[snp.seq[k]][snp.blocks[j][k].base]];

        // normalize score.
        double sum = 0.0;
        for(unsigned int j = 0; j < snp.probs.size(); j++)
            sum += snp.probs[j];
        if(sum == 0)
            error("sum = 0");
        for(unsigned int j = 0; j < snp.probs.size(); j++)
            snp.probs[j] /= sum;
    }
}

void do_snpcall(int argc, char **argv)
{
    prog_options options = parse_snpcall_options(argc, argv);

    cand_blocks_t cand_blocks;
    cand_vector_t snp_cands;

    // get the read length.
    unsigned int read_length = get_read_length(options.map_file.c_str());

	mapping_type t = get_mapping_type(options.map_file.c_str());

	if( t & COLOR )
	{
        // TODO: UGLY
        options.old_bin_cnt = 36;

		read_snp_cand_file(options.cand_file.c_str(), cand_blocks, snp_cands);

		array_type5 call_rates;
		boost::multi_array<int, 2> qual_bins;

        array_type5 call_counts;
        // If the call rates are not provided in a file, we need to load the background.
        if( options.call_file.size() == 0)
            fill_call_priors(options, call_counts);
        else
            read_call_counts_file(options.call_file.c_str(), call_counts);
        
		process_color_counts_to_rates(call_counts, call_rates, options.old_bin_cnt, options.new_bin_cnt, read_length, qual_bins);

		color_call_vector_t snp_calls;
		for(unsigned int i = 0; i < snp_cands.size(); i++)
		{
			SnpRecord &snp = snp_cands[i];
			snp_calls.push_back(ColorSnpCall(snp.id, snp.chr, snp.start, snp.end, snp.seq, snp.subs, snp.totl));
		}
		
		call_snps(options, cand_blocks, snp_calls, call_rates, read_length, qual_bins);

		boost::iostreams::filtering_ostream out;
		out.push(boost::iostreams::gzip_compressor());
		out.push(std::cout);
		boost::archive::binary_oarchive oa(out);
		oa << SNP_CALL_COL_FILE_TYPE;
		oa << SNP_CALL_COL_FILE_VERSION;
        oa << options.diploid;
		oa << snp_calls;
	}
	else // SEQUENCE
	{
        options.old_bin_cnt = 96;

		read_snp_cand_file(options.cand_file.c_str(), cand_blocks, snp_cands);

		array_type4 call_rates;
		boost::multi_array<int, 2> qual_bins;

        array_type4 call_counts;
        if( options.call_file.size() == 0 )
            fill_call_priors_seq(options, call_counts);
        else
            read_seq_call_counts_file(options.call_file.c_str(), call_counts);
        
    
		process_seq_counts_to_rates(call_counts, call_rates, options.old_bin_cnt, options.new_bin_cnt, read_length, qual_bins);

		seq_call_vector_t snp_calls;
		for(unsigned int i = 0; i < snp_cands.size(); i++)
		{
			SnpRecord &snp = snp_cands[i];
			snp_calls.push_back(SeqSnpCall(snp.id, snp.chr, snp.start, snp.end, snp.seq, snp.subs, snp.totl));
		}
		
		call_seq_snps(options, cand_blocks, snp_calls, call_rates, read_length, qual_bins);

		boost::iostreams::filtering_ostream out;
		out.push(boost::iostreams::gzip_compressor());
		out.push(std::cout);
		boost::archive::binary_oarchive oa(out);
		oa << SNP_CALL_SEQ_FILE_TYPE;
		oa << SNP_CALL_SEQ_FILE_VERSION;
        oa << options.diploid;
		oa << snp_calls;
	}

}

void call_snps(prog_options &options, cand_blocks_t &cand_blocks, color_call_vector_t &snp_cands, array_type5 &call_rates, unsigned int read_length, boost::multi_array<int, 2> &qual_bins)
{
    // For each snp candidate block generate list of model sets.
    // skip first dummy block.
    for(unsigned int i = 1; i < snp_cands.size(); i++)
    {
        ColorSnpCall &snp = snp_cands[i];

        // TODO: create aux. storage for coverage stats per snp site.
        if( options.diploid )
        {
            std::vector< std::vector<SnpModel> > elements;
            for(unsigned int j = 0; j < snp.end - snp.start; j++)
                elements.push_back(diploid_models[cswap[snp.seq[j]]]);
            snp.blocks = cartesian_product( elements );
        }
        else // haploid
            snp.blocks = cartesian_product( haploid_models, snp.end - snp.start );
        snp.probs.resize(snp.blocks.size(), 1.0);
    }

    ColorMappingParser parser(options.map_file.c_str());

    ColorRead read;
    while( parser.next_read(read) )
    {
        // ambig parse or unique parse.
        if( options.unique && read.num > 1 )
            continue;
        
        // find covered snps.
        // snp_map[ snp_id ] = ([read_id, cnt], ...) which mapping covers it.
        std::map<unsigned int, std::map<unsigned int, unsigned int> > snp_map;
        for(unsigned int j = 0; j < read.num; j++)
        {
            MapLoci &loc = read.map_locs[j];
            if(cand_blocks.find(loc.id) == cand_blocks.end())
            {
                if(options.ignore_missing) 
                   continue;
                else
                    error("The mapping references a sequence not in snp candidates. (Did you forget --ignore_missing?)");
            }
            for(unsigned int i = 0; i < read_length; i++)
                if( cand_blocks[loc.id][loc.pos + i] )
                    snp_map[cand_blocks[loc.id][loc.pos+i]][j]++;
        }
        
        // no need to consider reads that do not map to snp locations.
        if( !snp_map.size() )
            continue;

        // translate the quality score into our new score.
        for(unsigned int i = 0; i < read_length; i++)
            read.qual[i] = qual_bins[i][read.qual[i]];

        std::vector<double> align_scores;
        std::vector<std::string> dna_ref;
        
        align_scores.resize(read.num);
        dna_ref.resize(read.num);

        double sum_scores  = 0;
        double sum_scores2 = 0;

        // calculate all alignment scores to mapping locations using genome seq.
        for(unsigned int j = 0; j < read.num; j++)
        {
            MapLoci &loc = read.map_locs[j];

            dna_ref[j] = color_to_seq(loc.seq);

            double align_score = 1.0;

            for(unsigned int i = 0; i < dna_ref[j].size() - 1; i++)
                align_score *= call_rates[i+1][read.qual[i+1]][cswap[dna_ref[j][i]]][cswap[dna_ref[j][i+1]]][iswap[read.seq[i+1]]];

            if( align_score > 1 )
                error("align_score > 1");
            
            align_scores[j] = align_score;
            sum_scores  += align_score;
            sum_scores2 += align_score * align_score;
        }

        std::map<unsigned int, std::map<unsigned int, unsigned int> >::iterator mi;
        std::map<unsigned int, unsigned int>::iterator mi2;

        for(mi = snp_map.begin(); mi != snp_map.end(); ++mi)
        {
            // read covers snp #snp_id with its mapping #map_id.
            unsigned int snp_id = mi->first;

            // skip multiple mappings to same snp.
            if((mi->second).size() != 1)
                continue;

            mi2 = (mi->second).begin();

            // for blocks we will expect mi2 > 1.  do not filter on this.
            assert(mi2->second > 0);

            unsigned int map_id = mi2->first;
            
            MapLoci &m = read.map_locs[map_id];
            ColorSnpCall &snp = snp_cands[snp_id];

            unsigned int block_len = snp.end - snp.start;
            
            // make sure snp block must be spanned by read.
            // the snp block can not be at the 0th position, or overlap the ends.
            if ( m.pos >= snp.start || m.pos + m.seq.size() - 1 <= snp.end)
                continue;
            
            double ascore = align_scores[map_id];

            unsigned int idx;
            if( m.strand == '+' )
                idx = snp.start - m.pos;

            else // m.strand == '-' 
                idx = m.pos + read.seq.size() - snp.end;

            // mapped read does not span block.
            // TODO: should this ever be triggered?  Don't we catch all above?
            if(idx == 0 || idx + block_len > read.seq.size() - 2)
            {
                continue;
            }

            double new_score = ascore;
            // factor out align score due to block.
            for(unsigned int i = idx - 1; i < idx + block_len + 1; i++)
                new_score /= call_rates[i+1][read.qual[i+1]][cswap[dna_ref[map_id][i]]][cswap[dna_ref[map_id][i+1]]][iswap[read.seq[i+1]]];

            // blocks holds all of the permutations of models for this snp block.
            for(unsigned int i = 0; i < snp.blocks.size(); i++)
            {
                // notice: we have merged forward/reverse read handling.
                std::vector<SnpModel> models = snp.blocks[i];

                // NOTE: this makes a copy.  But we reverse it, so its good to not mess with the original.
                if( m.strand == '-' )
                {
                    models = snp.blocks[snp.blocks.size() - i - 1];
                    std::reverse(models.begin(), models.end());
                }

                double block_score = 1.0;
                    
                // TODO: should we have already filtered this?
                // TODO: should this include a + block_len?
                if( idx == 0 || idx > read.seq.size() - 2 )
                    continue;
                
                double s;

                // score first transition into block.
                s = 0.0;
                for(unsigned int k = 0; k < 4; k++)
                    s += call_rates[idx][read.qual[idx]][cswap[dna_ref[map_id][idx-1]]][ k ][iswap[read.seq[idx]]] * models[0].prob[k];
                
                block_score *= s;
                
                // score the transition leaving the block.
                unsigned int B = models.size();
                s = 0.0;
                for(unsigned int k = 0; k < 4; k++)
                    s += call_rates[idx+B][read.qual[idx+B]][ k ][cswap[dna_ref[map_id][idx+B]]][iswap[read.seq[idx+B]]] * models[B-1].prob[k];
                
                block_score *= s;

                // score intermediate range.
                for(unsigned int j = 1; j < models.size(); j++)
                {
                    s = 0.0;
                    for(unsigned int k = 0; k < 4; k++)
                    for(unsigned int l = 0; l < 4; l++)
                        s += call_rates[idx+j][read.qual[idx+j]][ k ][ l ][iswap[read.seq[idx+j]]] * models[j-1].prob[k] * models[j].prob[l];
                    if(s > 1)
                        error("s > 1");

                    block_score *= s;
                }

                if( block_score > 1.0 )
                    error("block score > 1.0");

                /*
                 * contribution = \sum_{map_loci} \frac{AS(read, loc | base)^2}{\sum_{map_loci} AS(read, loc | base)}
                 *
                 * The alignment scores squared [AS(...)^2] sum are accumulated into sum_scores2.
                 * The alignment score [AS(...)] sum is accumulated into sum_scores.
                 * To account for the probability given the base at this mapping loci, we need to take out the old alignment score for this location, 
                 * and add back the new score (new_score) which is the old score (ascore) with the block probability divided out, and factored back in (block_score).
                */
                double contribution = ( sum_scores2 - ascore * ascore + new_score * block_score * new_score * block_score ) / (sum_scores - ascore + new_score * block_score );

                if( contribution > 1.0 )
                    error("contribution > 1.0");

                snp.probs[i] *= contribution;

                if(snp.probs[i] > 1.0 )
                    error("probability for snp > 1.0");
            }

            snp.total_reads++;
            
            if( snp.total_reads % options.normalize_freq == 0 )
            {
                double sum = 0.0;
                for(unsigned int j = 0; j < snp.probs.size(); j++)
                    sum += snp.probs[j];
                for(unsigned int j = 0; j < snp.probs.size(); j++)
                    snp.probs[j] /= sum;
            }
        }
    }
}


int solve_col_snp(prog_options &options, ColorSnpCall &snp, std::vector<std::vector<double> > &M)
{
    if( snp.seq.empty() )
    {
        std::cerr << "WARNING: no read was found which spanned the entire snp block [" << snp.chr << " " << snp.start + 1 << " " << snp.end + 1 << "] of length " << snp.end - snp.start << std::endl;
        return -1;
    }
    

    double max_score = 0.0;
    int max_index = 0;

    for(unsigned int j = 0; j < snp.probs.size(); j++)
        if(snp.probs[j] > max_score)
        {
            max_score = snp.probs[j];
            max_index = j;
        }

    unsigned int block_len = snp.end - snp.start;

    // create a matrix to accumulate per base probs.
    M.resize(block_len);
    for(unsigned int j = 0; j < M.size(); j++)
        // TODO: use a models.size() variable or something here.
        if( !options.diploid )
        {
            M[j].resize(haploid_models.size());
            for(unsigned int k = 0; k < haploid_models.size(); k++)
                M[j][k] = 0.0;
        }
        else
        {
            M[j].resize(diploid_models[0].size());
            for(unsigned int k = 0; k < diploid_models[0].size(); k++)
                M[j][k] = 0.0;
        }

    // to get the indexes of the models from the snp_blocks we create a dummy integer matrix.
    std::vector<int> _dummy;
    if( !options.diploid )
    {
        for(unsigned int j = 0; j < haploid_models.size(); j++)
            _dummy.push_back(j);
    }
    else // diplod case.
    {
        for(unsigned int j = 0; j < diploid_models[0].size(); j++)
            _dummy.push_back(j);
    }


    std::vector<std::vector<int> > model_indexes = cartesian_product(_dummy, snp.end - snp.start);

    for(unsigned int j = 0; j < model_indexes.size(); j++)
        for(unsigned int k = 0; k < model_indexes[j].size(); k++)
            M[k][model_indexes[j][k]] += snp.probs[j];

    return max_index;
}

void print_col_snps(prog_options &options, color_call_vector_t &snp_cands)
{
    std::vector< std::vector<double> > M;

    // print numbers with decimals. 
    std::cout << std::setiosflags(std::ios::fixed) << std::setprecision(5);
    
    // For each snp candidate block generate list of model sets.
    for(unsigned int i = 1; i < snp_cands.size(); i++)
    {
        ColorSnpCall &snp = snp_cands[i];

        int max_index = solve_col_snp(options, snp, M);
        if( max_index < 0)
        {
            std::cerr << "WARNING: no read was found which spanned the entire snp block [" << snp.chr << " " << snp.start + 1 << " " << snp.end + 1 << "] of length " << snp.end - snp.start << std::endl;
            continue;
        }

        if(isnan(snp.probs[max_index]))
        {
            std::cerr << "WARNING: NaN encountered at " << snp.chr << " " <<  std::setw(10) << snp.start << " " << std::setw(10) << snp.end << " in probability calculation.  Try lowering normalization frequency (--freq | -f). Currently normalization frequency is " << options.normalize_freq << std::endl;

            if( options.only_show_changed )
                continue;

            for(unsigned int j = 0; j < snp.end - snp.start; j++)
            {
                std::cout << snp.chr << " " << std::setw(10) << snp.start + j + 1 << " " << snp.seq[j] << " " << snp.seq[j] << " : ";
                std::cout << snp.totl[j] << " : ";
                double epsilon = 0.001;
                for(unsigned int k = 0; k < M[j].size(); k++)
                {
                    if(!options.diploid)
                        std::cout << haploid_models[k].base << " " << M[j][k] << "\t";
                    else // options.diploid
                        std::cout << diploid_models[cswap[snp.seq[j]]][k].base << " " << M[j][k] << "\t";
                    if( M[j][k] < -epsilon )
                        error("Negative probability found for M value.");
                }
                std::cout << std::endl;
            }

            continue;
        }

        for(unsigned int j = 0; j < snp.end - snp.start; j++)
        {

            if( options.only_show_changed && snp.seq[j] == snp.blocks[max_index][j].base )
                continue;

            std::cout << snp.chr << " " << std::setw(12) << snp.start + j + 1 << " " << snp.seq[j] << " " << snp.blocks[max_index][j].base << " : ";
            
            std::cout << snp.totl[j] << " : ";
            
            double epsilon = 0.001;
            for(unsigned int k = 0; k < M[j].size(); k++)
            {
                if(!options.diploid)
                    std::cout << haploid_models[k].base << " " << M[j][k] << "\t";
                else // options.diploid
                    std::cout << diploid_models[cswap[snp.seq[j]]][k].base << " " << M[j][k] << "\t";
                if( M[j][k] < -epsilon )
                    error("Negative probability found for M value.");
            }
            std::cout << std::endl;
        }
    }
}

// TODO: modify all these programs to do SAM/PERM based on a templte parameter,
// since we have so much duplicated code.
void call_seq_snps(prog_options &options, cand_blocks_t &cand_blocks, seq_call_vector_t &snp_cands, array_type4 &call_rates, unsigned int read_length, boost::multi_array<int, 2> &qual_bins)
{
    // For each snp candidate block generate list of model sets.
    // skip first 'bad' block.
    for(unsigned int i = 1; i < snp_cands.size(); i++)
    {
        SeqSnpCall &snp = snp_cands[i];
        
        // TODO: create aux. storage for coverage stats per snp site.
        if( options.diploid )
            snp.blocks = diploid_models[cswap[snp.seq[0]]];
        else // haploid
            snp.blocks = haploid_models;

        snp.probs.resize(snp.blocks.size(), 1.0);
    }
    
    mapping_type t = get_mapping_type(options.map_file.c_str());

    if(t & PERM_FILE)
    {
        SequenceMappingParser parser(options.map_file.c_str());

        SeqRead read;
        while( parser.next_read(read) )
        {
            // ambig parse or unique parse.
            if( options.unique && read.num > 1 )
                continue;
            
            // find covered snps.
            // snp_map[ snp_id ] = ([read_id, cnt], ...) which mapping covers it.
            std::map<unsigned int, std::map<unsigned int, unsigned int> > snp_map;
            for(unsigned int j = 0; j < read.num; j++)
            {
                MapLoci &loc = read.map_locs[j];
                if(cand_blocks.find(loc.id) == cand_blocks.end())
                {
                    if(options.ignore_missing) 
                        continue;
                    else
                        error("The mapping references a sequence not in snp candidates. (Did you forget --ignore_missing?)");
                }
                for(unsigned int i = 0; i < read_length; i++)
                    if( cand_blocks[loc.id][loc.pos + i] )
                        snp_map[cand_blocks[loc.id][loc.pos+i]][j]++;
            }
            
            // no need to consider reads that do not map to snp locations.
            if( !snp_map.size() )
                continue;

            // translate the quality score into our new score.
            for(unsigned int i = 0; i < read_length; i++)
                read.qual[i] = qual_bins[i][read.qual[i]];

            std::vector<double> align_scores;
            
            align_scores.resize(read.num);

            double sum_scores  = 0;
            double sum_scores2 = 0;

            // calculate all alignment scores to mapping locations using genome seq.
            for(unsigned int j = 0; j < read.num; j++)
            {
                MapLoci &loc = read.map_locs[j];

                double align_score = 1.0;

                for(unsigned int i = 0; i < loc.seq.size(); i++)
                    align_score *= call_rates[i][read.qual[i]][cswap[loc.seq[i]]][cswap[read.seq[i]]];

                if( align_score > 1 )
                    error("align_score > 1");
                
                align_scores[j] = align_score;
                sum_scores  += align_score;
                sum_scores2 += align_score * align_score;
            }

            std::map<unsigned int, std::map<unsigned int, unsigned int> >::iterator mi;
            std::map<unsigned int, unsigned int>::iterator mi2;

            for(mi = snp_map.begin(); mi != snp_map.end(); ++mi)
            {
                // read covers snp #snp_id with its mapping #map_id.
                unsigned int snp_id = mi->first;

                // skip multiple mappings to same snp.
                if((mi->second).size() != 1)
                    continue;

                mi2 = (mi->second).begin();

                // TODO ???
                // for blocks we will expect mi2 > 1.  do not filter on this.
                assert(mi2->second > 0);

                unsigned int map_id = mi2->first;
                
                MapLoci    &m   = read.map_locs[map_id];
                SeqSnpCall &snp = snp_cands[snp_id];

                double ascore = align_scores[map_id];

                unsigned int idx;
                if( m.strand == '+' )
                    idx = snp.start - m.pos;
                else // m.strand == '-' 
                    idx = m.pos + read.seq.size() - snp.end;

                // factor out align score due to block.
                double new_score = ascore / call_rates[idx][read.qual[idx]][cswap[m.seq[idx]]][cswap[read.seq[idx]]];

                // blocks holds all of the permutations of models for this snp block.
                for(unsigned int i = 0; i < snp.blocks.size(); i++)
                {
                    // notice: we have merged forward/reverse read handling.
                    SnpModel model = snp.blocks[i];

                    // NOTE: this makes a copy.
                    // we do not revese since its a single base.
                    if( m.strand == '-' )
                        model = snp.blocks[snp.blocks.size() - i - 1];

                    double block_score = 0.0;
                    // score intermediate range.
                    for(unsigned int k = 0; k < 4; k++)
                        block_score += call_rates[idx][read.qual[idx]][ k ][cswap[read.seq[idx]]] * model.prob[k];
                    if(block_score > 1)
                        error("block_score > 1");
                    /*
                     * contribution = \sum_{map_loci} \frac{AS(read, loc | base)^2}{\sum_{map_loci} AS(read, loc | base)}
                     *
                     * The alignment scores squared [AS(...)^2] sum are accumulated into sum_scores2.
                     * The alignment score [AS(...)] sum is accumulated into sum_scores.
                     * To account for the probability given the base at this mapping loci, we need to take out the old alignment score for this location, 
                     * and add back the new score (new_score) which is the old score (ascore) with the block probability divided out, and factored back in (block_score).
                    */
                    double contribution = ( sum_scores2 - ascore * ascore + new_score * block_score * new_score * block_score ) / (sum_scores - ascore + new_score);

                    if( contribution > 1.0 )
                        error("contribution > 1.0");

                    snp.probs[i] *= contribution;

                    if(snp.probs[i] > 1.0 )
                        error("probability for snp > 1.0");
                }

                snp.total_reads++;
                
                if( snp.total_reads % options.normalize_freq == 0 )
                {
                    double sum = 0.0;
                    for(unsigned int j = 0; j < snp.probs.size(); j++)
                        sum += snp.probs[j];
                    for(unsigned int j = 0; j < snp.probs.size(); j++)
                        snp.probs[j] /= sum;
                }


            }
        }
    }
    else if(t & SAM_FILE)
    {
        std::map<std::string, std::string> seqs;
        parse_fasta(options.fasta_file.c_str(), seqs);
        SamSequenceMappingParser parser(options.map_file.c_str(), seqs, t);

        SeqRead read;
        while( parser.next_read(read) )
        {
            // ambig parse or unique parse.
            if( options.unique && read.num > 1 )
                continue;
            
            // find covered snps.
            // snp_map[ snp_id ] = ([read_id, cnt], ...) which mapping covers it.
            std::map<unsigned int, std::map<unsigned int, unsigned int> > snp_map;
            for(unsigned int j = 0; j < read.num; j++)
            {
                MapLoci &loc = read.map_locs[j];
                if(cand_blocks.find(loc.id) == cand_blocks.end())
                {
                    if(options.ignore_missing) 
                        continue;
                    else
                        error("The mapping references a sequence not in snp candidates. (Did you forget --ignore_missing?)");
                }
                for(unsigned int i = 0; i < read_length; i++)
                    if( cand_blocks[loc.id][loc.pos + i] )
                        snp_map[cand_blocks[loc.id][loc.pos+i]][j]++;
            }
            
            // no need to consider reads that do not map to snp locations.
            if( !snp_map.size() )
                continue;

            // translate the quality score into our new score.
            for(unsigned int i = 0; i < read_length; i++)
                read.qual[i] = qual_bins[i][read.qual[i]];

            std::vector<double> align_scores;
            
            align_scores.resize(read.num);

            double sum_scores  = 0;
            double sum_scores2 = 0;

            // calculate all alignment scores to mapping locations using genome seq.
            for(unsigned int j = 0; j < read.num; j++)
            {
                MapLoci &loc = read.map_locs[j];

                double align_score = 1.0;

                for(unsigned int i = 0; i < loc.seq.size(); i++)
                    align_score *= call_rates[i][read.qual[i]][cswap[loc.seq[i]]][cswap[read.seq[i]]];

                if( align_score > 1 )
                    error("align_score > 1");
                
                align_scores[j] = align_score;
                sum_scores  += align_score;
                sum_scores2 += align_score * align_score;
            }

            std::map<unsigned int, std::map<unsigned int, unsigned int> >::iterator mi;
            std::map<unsigned int, unsigned int>::iterator mi2;

            for(mi = snp_map.begin(); mi != snp_map.end(); ++mi)
            {
                // read covers snp #snp_id with its mapping #map_id.
                unsigned int snp_id = mi->first;

                // skip multiple mappings to same snp.
                if((mi->second).size() != 1)
                    continue;

                mi2 = (mi->second).begin();

                // TODO ???
                // for blocks we will expect mi2 > 1.  do not filter on this.
                assert(mi2->second > 0);

                unsigned int map_id = mi2->first;
                
                MapLoci    &m   = read.map_locs[map_id];
                SeqSnpCall &snp = snp_cands[snp_id];

                double ascore = align_scores[map_id];

                unsigned int idx;
                if( m.strand == '+' )
                    idx = snp.start - m.pos;
                else // m.strand == '-' 
                    idx = m.pos + read.seq.size() - snp.end;

                // factor out align score due to block.
                double new_score = ascore / call_rates[idx][read.qual[idx]][cswap[m.seq[idx]]][cswap[read.seq[idx]]];

                // blocks holds all of the permutations of models for this snp block.
                for(unsigned int i = 0; i < snp.blocks.size(); i++)
                {
                    // notice: we have merged forward/reverse read handling.
                    SnpModel model = snp.blocks[i];

                    // NOTE: this makes a copy.
                    // we do not revese since its a single base.
                    if( m.strand == '-' )
                        model = snp.blocks[snp.blocks.size() - i - 1];

                    double block_score = 0.0;
                    // score intermediate range.
                    for(unsigned int k = 0; k < 4; k++)
                        block_score += call_rates[idx][read.qual[idx]][ k ][cswap[read.seq[idx]]] * model.prob[k];
                    if(block_score > 1)
                        error("block_score > 1");
                    /*
                     * contribution = \sum_{map_loci} \frac{AS(read, loc | base)^2}{\sum_{map_loci} AS(read, loc | base)}
                     *
                     * The alignment scores squared [AS(...)^2] sum are accumulated into sum_scores2.
                     * The alignment score [AS(...)] sum is accumulated into sum_scores.
                     * To account for the probability given the base at this mapping loci, we need to take out the old alignment score for this location, 
                     * and add back the new score (new_score) which is the old score (ascore) with the block probability divided out, and factored back in (block_score).
                    */
                    double contribution = ( sum_scores2 - ascore * ascore + new_score * block_score * new_score * block_score ) / (sum_scores - ascore + new_score);

                    if( contribution > 1.0 )
                        error("contribution > 1.0");

                    snp.probs[i] *= contribution;

                    if(snp.probs[i] > 1.0 )
                        error("probability for snp > 1.0");
                }

                snp.total_reads++;
                
                if( snp.total_reads % options.normalize_freq == 0 )
                {
                    double sum = 0.0;
                    for(unsigned int j = 0; j < snp.probs.size(); j++)
                        sum += snp.probs[j];
                    for(unsigned int j = 0; j < snp.probs.size(); j++)
                        snp.probs[j] /= sum;
                }
            }
        }
    }
    else
        error("Bad File type");

}

void print_seq_snps(prog_options &options, seq_call_vector_t &snp_cands)
{
    std::cout << std::setiosflags(std::ios::fixed) << std::setprecision(5);

    multiply_seq_snp_calls_by_priors(options, snp_cands);

    // For each snp candidate block generate list of model sets.
    for(unsigned int i = 1; i < snp_cands.size(); i++)
    {
        SeqSnpCall &snp = snp_cands[i];

        if( snp.seq.empty() )
        {
            std::cerr << "WARNING: no read was found which spanned the snp [" << snp.chr << " " << snp.start + 1 << "]" << snp.end - snp.start << std::endl;
            continue;
        }

        int max_index = std::max_element(snp.probs.begin(), snp.probs.end()) - snp.probs.begin();

        if( options.only_show_changed && snp.seq[0] == snp.blocks[max_index].base )
            continue;

        std::cout << snp.chr << " " << std::setw(10) << snp.start + 1 << " " << snp.seq << " " << snp.blocks[max_index].base << " : ";

        std::cout << snp.totl[0] << " : ";
        
        if(!options.diploid)
        {
            for(unsigned int j = 0; j < haploid_models.size(); j++)
                std::cout << haploid_models[j].base << " " << snp.probs[j] << "\t";
        }
        else
        {
            for(unsigned int j = 0; j < diploid_models[cswap[snp.seq[0]]].size(); j++)
                std::cout << diploid_models[cswap[snp.seq[0]]][j].base << " " << snp.probs[j] << "\t";
        }
        std::cout << std::endl;
    }
}
